function F_star = DCCS_rcl(df_sim,coefs_sim,dfull_sim, nmarket)

brndfe = df_sim.brndid == unique(df_sim.brndid)';
montfe = df_sim.montid == unique(df_sim.montid)';
prodfe = df_sim.prodid == unique(df_sim.prodid)';

% Random sample X of markets
sampled_cdid = (1:nmarket)';
% sampled_cdid = randsample(unique_cdid, nmarket);

sel = ismember(df_sim.cdid_sim, sampled_cdid);
df_sim = df_sim(sel, :);
dfull_sim = dfull_sim(sel, :);

% data parameters
fesd     = [prodfe(sel, :) montfe(sel, :)];
x        = [df_sim.p_jt fesd];
xnonlin  = [df_sim.p_jt ones(size(df_sim.p_jt)) df_sim.calor];
xi       = [df_sim.xi];
theta    = [coefs_sim.alpha; coefs_sim.sigma_D];
Pi       = coefs_sim.Pi;
rho      = coefs_sim.rho;

% model parameters
fes_M       = [brndfe(sel, :) montfe(sel,1:end-1)];
x_M         = [df_sim.p_jt fes_M];
theta_M0 = [coefs_sim.alpha; Pi(2); Pi(3)];
theta_M = theta_M0;

% Make wide by market
x_mrkwide = zeros(size(prodfe, 2), size(x, 2), nmarket);
xi_mrkwide = zeros(size(prodfe, 2), size(xi, 2), nmarket);
xnonlin_mrkwide = zeros(size(prodfe, 2), size(xnonlin, 2), nmarket);
dfull_sim_mrkwide = zeros(size(prodfe, 2), size(dfull_sim, 2), nmarket);
in_samp_mrkwide = zeros(size(prodfe, 2), size(dfull_sim, 2), nmarket);

for m = 1:nmarket
    sel = (df_sim.cdid_sim == m);
    x_mrkwide(1:size(x(sel,:),1),1:size(x(sel,:),2),m) = x(sel,:);
    xi_mrkwide(1:size(xi(sel,:),1),1:size(xi(sel,:),2),m) = xi(sel,:);
    xnonlin_mrkwide(1:size(xnonlin(sel,:),1),1:size(xnonlin(sel,:),2),m) = xnonlin(sel,:);
    dfull_sim_mrkwide(1:size(dfull_sim(sel,:),1),1:size(dfull_sim(sel,:),2),m) = dfull_sim(sel,:);
    in_samp_mrkwide(1:size(dfull_sim(sel,:),1),1:size(dfull_sim(sel,:),2),m) = 1;
end

% Calculate data shares for data
delta   = pagemtimes(x_mrkwide, reshape(theta, [], 1, 1)) + xi_mrkwide;
[share_main, ~, der] = rcnl_der_price_mktwide( ...
    theta(1),           ...
    delta,              ...
    Pi,                 ...
    rho,                ...
    xnonlin_mrkwide,     ...
    dfull_sim_mrkwide,   ...
    in_samp_mrkwide      ...
);

f_obj = @(theta_val) ObjFun(theta_val,xnonlin_mrkwide,dfull_sim_mrkwide, in_samp_mrkwide, der, share_main, nmarket);

   options = optimoptions('fminunc',     ...
        'Display', 'iter-detailed',       ...
        'MaxIterations',1000,              ...
        'Algorithm', 'quasi-newton',    ...
        'StepTolerance', 1e-10,           ...
        'SpecifyObjectiveGradient', true, ...
        'OptimalityTolerance',1e-6,       ...
        'checkGradients',false            ...
    );

% Perform optimization
[theta_star] = fminunc(f_obj, theta_M0, options);

% Return optimized value
[F_star, dF_star] = ObjFun(theta_star,xnonlin_mrkwide,dfull_sim_mrkwide, in_samp_mrkwide, der, share_main, nmarket);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [F,dF] = ObjFun(theta_M,xnonlin_mrkwide,dfull_sim_mrkwide, in_samp_mrkwide, der, share_main, nmarket)

Pi_M = [0; theta_M(end-1:end)];
rho_M = 0;
fe_list = theta_M(1);

% Calculate delta for model
[Delta_RCL, dDelta_RCL] = rcnl_delta( ...
    share_main,         ...
    Pi_M,            ...
    rho_M,           ...
    xnonlin_mrkwide,          ...
    dfull_sim_mrkwide          ...
);

[~, share_main_ijt_M, der_M] = rcnl_der_price_mktwide( ...
    fe_list(1),           ...
    Delta_RCL,              ...
    Pi_M,                 ...
    rho_M,                ...
    xnonlin_mrkwide,     ...
    dfull_sim_mrkwide,   ...
    in_samp_mrkwide      ...
);

% gradient intialize
grad = zeros(size(share_main_ijt_M,1), size(share_main_ijt_M,1), size(share_main_ijt_M,3), 2);

for i = 1:size(theta_M,1)
    % for linear parameters
    if i == 1
        if fe_list(1) ~= 0
            grad(:,:,:,i) = (1/fe_list(1)).*der_M;
        else 
            grad(:,:,:,i) = ones(size(grad(:,:,:,i)));
        end
    % for non-linear parameters
    else
        [~,ds_dtheta] = f_calc_gradient_mktwide( ...
            xnonlin_mrkwide(:,i,:), ...
            dfull_sim_mrkwide,   ...
            share_main_ijt_M, ...
            dDelta_RCL(:,i,:) ...
        );
        % ds_dtheta = repmat(ds_dtheta, [1, 500, 1]);

        [J, ns, ~]    = size(share_main_ijt_M);
    
        eyeJ = eye(J);
        eye3D = repmat(eyeJ, 1, 1, nmarket);
        
        own   = sum(fe_list(1).*(1 - share_main_ijt_M) .* ds_dtheta, 2);
    
        lhs = -fe_list(1)*pagemtimes(share_main_ijt_M, pagetranspose(ds_dtheta));
        lhs(logical(eye3D)) = own;
        rhs = fe_list(1)*pagemtimes(ds_dtheta, pagetranspose(share_main_ijt_M));
    
        grad(:,:,:,i) = 1/ns .* (lhs - rhs);
    end    
end

F = 1/nmarket*sum(sum(sum((der-der_M).^2)));
dF = zeros(length(theta_M),1);
for p = 1:length(theta_M)
dF(p) = -1/nmarket*sum(sum(sum(2.*(der-der_M).*grad(:,:,:,p))));
end

end
